/**
 * The binomial coefficient C(10^18,10^9) is a number with more than 9 billion
 * digits.
 *
 * Let M(n,k,m) denote the binomial coefficient C(n,k) modulo m.
 *
 * Calculate M(10^18,10^9,p*q*r) for 1000 < p < q < r < 5000 and p,q,r prime.
 */

#include <iostream>
#include <algorithm>
#include "euler/binomial.hpp"
#include "euler/prime_table.hpp"
#include "euler/gcd.hpp"
#include "euler.h"

BEGIN_PROBLEM(365, solve_problem_365)
	PROBLEM_TITLE("A huge binomial coefficient")
	PROBLEM_ANSWER("162619462356610313")
	PROBLEM_DIFFICULTY(1)
	PROBLEM_FUN_LEVEL(1)
	PROBLEM_TIME_COMPLEXITY("n^3")
	PROBLEM_SPACE_COMPLEXITY("n^2")
	PROBLEM_KEYWORDS("binomial")
END_PROBLEM()

#if 0
// Todo: optimize this routine to avoid modular multiplication.
static long long solve_congruence_3(int p1, int p2, int p3, int r1, int r2, int r3)
{
	int p[3] = { p1, p2, p3 };
	int r[3] = { r1, r2, r3 };

	long long m = (long long)p1 * p2 * p3;
	euler::residue<long long> x(0, m);
	for (int i = 0; i < 3; i++)
	{
		if (r[i] != 0)
		{
			x += r[i] * (m / p[i] * euler::modinv((int)(m/p[i]%p[i]), p[i]));
		}
	}
	return x.value();
}
#endif

static void solve_problem_365()
{
	const long long N = 1000000000000000000LL;
	const long long K = 1000000000LL;

	// List all primes between 1000 and 5000.
	euler::prime_table<int> ptable(5000);
	std::vector<int> P(ptable.lower_bound(1000), ptable.end());

	// Compute the residue of C(10^18, 10^9) modulo each prime.
	std::vector<int> R(P.size());
	for (size_t i = 0; i < P.size(); i++)
	{
		R[i] = euler::modbinom(N, K, P[i]);
	}

	// For each pair of primes (p,q), compute q*modinv(q,p) and p*modinv(p,q).
	int count = (int)P.size();
	std::vector<int> cache(count * count);
#if 0
	for (int i = 0; i < count; i++)
	{
		for (int j = 0; j < count; j++)
		{
			if (i != j && R[i] != 0)
			{
				cache[i*count+j] = P[j]*euler::modinv(P[j],P[i]);
			}
		}
	}
#else
	for (int i = 0; i < count; i++)
	{
		for (int j = i + 1; j < count; j++)
		{
			if (R[i] != 0 || R[j] != 0)
			{
				int x = 0, y = 0;
				euler::egcd(P[i], P[j], x, y);
				if (R[i])
					cache[i*count+j] = P[j] * (y < 0? y + P[i] : y);
				if (R[j])
					cache[j*count+i] = P[i] * (x < 0? x + P[j] : x);
			}
		}
	}
#endif

	// For each prime triple p1 < p2 < p3, compute C(N,K) mod p1*p2*p3
	// using the Chinese remainder theorem.
	long long result = 0;
	for (int i = 0; i < count; i++)
	{
		for (int j = i + 1; j < count; j++)
		{
			for (int k = j + 1; k < count; k++)
			{
#if 1
				long long m = (long long)P[i]*P[j]*P[k];
				long long t = 0;
				if (R[i])
					t += R[i] * ((long long)cache[i*count+j] * cache[i*count+k] % m);
				if (R[j])
					t += R[j] * ((long long)cache[j*count+k] * cache[j*count+i] % m);
				if (R[k])
					t += R[k] * ((long long)cache[k*count+i] * cache[k*count+j] % m);
				if (t)
					result += t % m;
#else
				result += solve_congruence_3(P[i],P[j],P[k],R[i],R[j],R[k]);
#endif
			}
		}
	}
	std::cout << result << std::endl;
}
